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We investigate the possibility of reconstructing the initial spectrum of density fluctu- 
ations from the cosmic microwave background (CMB) anisotropy. As a first step toward 
this program, we consider a spatially flat, CDM dominated universe. In this case, it is 
shown that, with a good accuracy, the initial spectrum satisfies a first order differen- 
tial equation with the source determined by the CMB angular correlation function. The 
equation is found to contain singularities arising from zeros of the acoustic oscillations 
in the transfer functions. Nevertheless, we find these singularities are not fatal, and the 
equation can be solved nicely. We test our method by considering simple analytic forms 
for the transfer functions. We find the initial spectrum is reproduced within 5% accuracy 
even for a spectrum that has a sharp spike. 

PACS: 95.30.-k; 98.80.Es 



I. INTRODUCTION 

Determination of the primordial spectrum of the density fluctuations is one of the most important 
issues in modern cosmology. The cosmic microwave background (CMB) anisotropy provides us with a 
great deal of information of the primordial fluctuations, and it is considered to be a powerful tool for 
studying the early universe |l| . This is because the physical processes that determine the CMB anisotropy 
are described by linear perturbation theory, and they have been well understood 

Although the recent anisotropy observations are consistent with a flat universe with a scale-invariant 
initial power spectrum, this does not exclude the possibilities of other models. In particular, we do not 
know how much the shape of the initial power spectrum is constrained. In most of previous investigations, 
when cosmological model parameters are estimated from the observational data by likelihood analysis, 
the initial spectrum is assumed to have a power-law shape M. It is true that a conventional slow- 
roll inflation model that has now become a 'standard model', gives a power-law spectrum which 
is almost scale-invariant However, when analyzing the observed CMB anisotropy, it is much more 
desirable, and probably much healthier, to constrain the initial spectrum solely by observed data without 
any theoretical prejudices. For example, even within the context of inflationary cosmology, a variety of 
generation mechanisms for non-scale-invariant perturbations have been proposed [Q. In this connection, 
recently several authors have discussed extraction of non-power law features from the CMB observations 
1^, where the initial spectrum is allowed to have a piece- wise power-law shape. 

In this paper, we approach this issue in an entirely different way, namely, by formulating an inverse 
problem as faithful as possible. Such an approach will be eventually needed if we seriously want to 
constrain the initial power spectrum solely from observations of the CMB anisotropy. An approach to 
this inversion problem has been discussed recently 

As a first step, we consider a simple situation in which the transfer functions that relate the input 
power spectrum P{k) of the gravitational potential to the output CMB angular correlation function C{9) 
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are given analytically. This is certainly a toy model. However, it has almost all the essential features a 
realistic model would have. In particular, unlike |9), our model takes account of not only the Sachs- Wolfe 
(SW) effect [nO[ but also the Doppler effect . The latter, which gives rise to zero points in the transfer 
functions, is tEe main cause of the difficulty in this inversion problem. 

The advantage of adopting this simple situation is that our method of inversion which we shall develop 
below may be easily tested at various stages of calculations. Since our primary concern here is to formulate 
the inversion problem, we fix the cosmological parameters and do not study the dependence of P{k) upon 
them. 

The paper is organized as follows. In Sec. II, we give the basic equations that relate the primordial 
spectrum P{k) with the angular correlation function C{9). In Sec. Ill, under some reasonable assump- 
tions, we derive a differential equation for P{k) and develop a method to solve it. Then we test our 
method by applying it to several spectral shapes. We find our method is applicable even in the case of a 
spectrum with a sharp spike. 



II. BASIC EQUATIONS 

The angular correlation function of the CMB, C{9) is defined as 

= (8(71)6(72)), COS0-71-72 (1) 

where 6(7) — {AT/T){-f) is the temperature fluctuation in the direction 7 and the average is taken over 
all angles and all spatial positions. C{d) is expanded in the Legendre polynomials and related to the 
angular power spectrum, Ci, as 

2/ 4- 1 

^(^)=E^^'^'(c°«^)- (2) 
I 

C{9) or Ci are the fundamental observables of the CMB anisotropy. Our purpose is to reconstruct the 
initial power spectrum from the above observables. We calculate in the Newtonian gauge and follow the 
notation of 

Each Fourier mode of temperature perturbations, 6(?7, fc), obeys the Boltzmann equation JTl[ |, 

e + ikfi{e + *) = -$ + f[eo - e - ^e2P2{f^) - ^^J-Vbl (3) 

where the dot denotes a derivative with respect to the conformal time 77, k is the comoving wave number, 
fj, = k~^k • 7, f is the differential Thomson optical depth and Vb is the bulk velocity of baryons. "if and $ 
are the gauge-invariant Newtonian potential and spatial curvature perturbation, respectively [ pl] |. Here, 
we decompose Q{r], k, /i) into the multipole moments, 

&{7^,k,fl) = J2i-^yQliv,k)Pli^^), (4) 
I 

where Qiit], k) is the l-th multipole moment of Q{rj, k, /i). Integrating Eq. (^), we obtain 

(9 + ^){vo,k, m) - ^ {[Oo + * - z/il4]V(77) + (* - $)£-"('') } e''=^("-"«)dr/, (5) 

where rjo is the conformal time today and V{ri) is the visibility function given by 

rno 

Viv) = r(,7)e--('') ; ri^) = / T{i^')drj'. (6) 



J J] 

We have neglected the quadrupole term on the right hand side of Eq. (||) since its contribution is negligible 
in the tight coupling approximation. The visibility function has a sharp peak around the last scattering 
time 77* so that we assume that recombination occurs instantaneously at 77 = ry*. Then, the multipole 
moments of each k mode is approximately given by 

9,(770, A:) = (60 + ^){7j,,k){2l + l)n{kd) + 61(77,, fc)(2/ + l)j,'(fcd) 

+ {21 + 1)1 dr;' — (*(77',fc)-$(77',fc))jKfc%- W), (7) 



where d = jiq — rj^ is a conformal distance from the present epoch to the last scattering surface (LSS). 
Equation (M) shows that there are four sources for the temperature anisotropy, namely, the intrinsic 
temperature variation, the Sachs- Wolfe (SW) effect which is caused by the static gravitational potential 
at the last scattering surface and is dominant on large scales, the Doppler effect due to the fluid bulk 
velocity, and the integrated Sachs- Wolfe (ISW) effect due to the time variation of the metric perturbation 
after decoupling. If the recombination occurs when the universe is matter-dominated, ^ and $ remain 
constant in time after decoupling until a stage at which either the curvature term or the cosmological 
constant term becomes significant. In particular, in the Einstein-de Sitter universe, they remain constant 
until today, hence the ISW effect is absent. 
Conventionally, C; is expressed as 
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From Eqs. M) and (H), we find 



Ci 
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dk k^ 



{Qo + ^)(Tl*)ji(kd) + QMj'i{kd) 



(9) 



The angular correlation function is calculated from Eq. (0). Here, we define r by 



Q 

r = 2c? sin — . 
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(10) 



This is the spatial distance between two points on the last scattering surface which are observed with 
the angular separation 6. Since the thickness of the LSS is neglected in Eq. (^, there is a one-to-one 
correspondence between the observed temperature anisotropy and the perturbation variables on the LSS. 
Using the relation, 



Y,i2l + l)Pi{cos9)jfikd) = 
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the angular correlation function C(r) is given by 
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In the linear perturbation theory, different k modes do not couple with each other but evolve indepen- 
dently. The initial condition for each k mode is directly reflected in the perturbation variables through 
transfer functions. The variables, Qq, ©i and 4", in the integrand of Eq. (O) are therefore linearly related 

to the initial condition and we can generally express the relation between C(r) and the initial spectrum 
as 

/"OO 

C{r) = / K{k,r)P{k)dk, (13) 
Jo 

where we normalize the initial condition in terms of the curvature perturbation $(r/ — 0,k) and define 
P{k) = (|$(0,A:)p). In this case, introducing the transfer functions /(fc) and ^(fc) defined by 
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(eo + *)(r7*,fc) = /(^*,A:)$(0,fc), 



K{k,r) can be written as 



(14) 



K{k,r) = 



where ai{k) are given by 
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Our purpose is to solve the integral equation ( [l3) f or P{k) with the kernel given by Eq. ([l5|). When 
there is only the monopole term, i.e., Qi — 0, Eq. (1 



Cir) = 



27r2 



131) takes the familiar form, 
' dkk^ /2(fc)F(fc)-''''^'' 
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In this case, using the Fourier sine formula, we obtain 

An 



P{k) 



P{k)k Jo 



dr r C{r) sin kr. 



(17) 



(18) 



It must be noted that P{k) has singularities on small scales since f{k) and g{k) are oscillatory functions 
reflecting the acoustic oscillations of the density perturbations. These singularities are inevitable as long 
as we take this approach and the one-to-one correspondence between C{r) and P{k) breaks down at the 
singularities. However, we shall find in the next section that there is a way to resolve this difficulty. 



III. INVERSION METHOD 

In the above discussion, we have tacitly assumed that r runs from zero to infinity. In reality, however, 
r is bounded in the finite range < r < 2d. Furthermore, it is observationally impossible to determine 
C(r) on large scales due to the statistical ambiguity, i.e., the cosmic variance. However, the scales which 
we are interested in are r <^ d and it is expected that modes with k > 1/d have little effect on these 
scales. We therefore neglect the terms proportional to 1/d and l/d^ in Eq. (|l^). In this limit, we have 

C{r)^-^( dk^^ {F{k)k^r^sm{kr)-G{k)kr cos kr + G{k)suikr} , (19) 
27r Jo kr 

where we have replaced /^(fc) and g'^{k) with F{k) and G{k), respectively, for notational simplicity. 
Multiplying by r^ and integrating by parts twice, we may re-express Eq. ( [l9| ) as 

^'^W^^y^ dk ^-^{kF{k)P{k)) + —{G{k)P{k)) + -Gik)P{k)jsinkr, (20) 

where we have assumed that the boundary terms vanish and the integral converges. Employing the 
Fourier sine formula, we obtain 



A 



kF{k)P"{k)+{-2F - 2kF' + G)P'{k) 



-(C - 2F' - kF" + G/k)P{k) = 47r / 



r^Cir)smkrdr. (21) 



This is a second-order ordinary differential equation in whicli the source term carries the information 
of C(r). It is of course necessary to fix the boundary conditions for P{k) in order to solve Eq. (^l|). 
However, we have no means to determine these conditions. The only knowledge we have is that The 
boundary conditions at A: = and k = oo are restricted from the convergence of the integral in Eq. (p^. 

In order to avoid this difficulty, we consider the reduction of the order of the differential equation. If 
we take some appropriate combination of C(r) and its derivative, we can derive a differential equation of 
lower order. The simplest combination is given by 

1 r°° 

C{r) = 3rC{r) + r'^C'{r) = — dk P(k) {F{k)k^r cos kr + (2F{k) + G(k))k sin kr} . (22) 





Note that the power of r in the integrand is reduced. No linear combination can reduce the order of the 
differential equation lower than C{r). Integrating by parts, we obtain 

/>oo 

- Fk^P' + {-F'k + G)kP ^ in C{r)smkrdr. (23) 

Jo 

This is a first order differential equation. Wc discuss a method for solving this equation in the rest of 
this section. r_. 

We now give the expressions for f{k) and g{k) explicitly in order to examine the properties of Eq. (|23|), 
and to verify if thus obtained solution correctly reproduces the original spectrum. For a given model, 
f{k) and g{k) are determined by the coupled Einstein-fluid equations which are to be solved numerically. 
However, in order to understand the property of Eq. ( p3| ) and to establish a method which we can apply 
to general cases, we start with a toy model in which jjk) and g{k) are given analytically. 

Before recombination, photons are tightly coupled with baryons through the Thomson scattering. With 
the tight coupling approximation, we can expand the Boltzmann equation for photons together with the 
continuity and Euler equations for baryons in the Thomson scattering time To the leading order, the 
dynamics of the photon-baryon fluid is described by the following simple equation for each Fourier mode 

al + K 

where Cs is the sound speed of the photon-baryon fluid, R = ipb/ip^, and T{r]) is the driving term due 
to the gravitational potential, 

Hv)-~^--T^^-^^- (25) 
al + R 3 

The dipole moment is obtained from the continuity equation, 

fcOi = -3(6)0 + *). (26) 

Note that Eqs. (|3) and (^) are applicable not only to the Einstein-de Sitter universe but also to other 
cosmological models. 

Here we solve Eqs. (|4|) and ( p6| ) under the following assumptions: 

(i) The gravitational potential fluctuation ^'(r;) and the spatial curvature fluctuation ^(rj) are always 
independent of time, i.e., ^I* = $ = 0, and ^I* = — $. 

(ii) The baryon denisty is negligible, i.e., R — 0. 

(iii) The perturbation is adiabatic. 

The assumption (i) is valid only on large and intermediate scales, but fails on small scales (i.e., on scales 
much smaller than the Hubble horizon scale at decoupling) since the acoustic oscillations tend to destroy 
the gravitational potential on such scales. The assumption (ii) is, of course, quite unrealistic, while (iii) 
is a natural consequence of inflationary cosmology. However, these assumptions do not alter the essential 



features of the CMB anisotropy spectrum. The difficulty of solving Eq. (|2^) is caused by the acoustic 
oscillations of temperature fluctuations which give rise to zero points in the transfer functions / and g, 
and this is a common feature of the CMB anisotropy regardless of cosmological models. Furthermore, we 
do not take account of the diffusion damping which becomes significant on small scales, since this 
effect do not significantly change the structure of Eq. (|2^) either. 

We solve Eqs. (|2j) and ( p6| ) in this toy model. The solutions are given by 

[60 + *] (77) - i$(0)cos(fcc,77), 

81(77) =c,$(0)sin(fcc,77), (27) 

where Cs — l/\/3 and we take the adiabatic initial condition as Oo(0) = ^'(0)/3 = — $(0)/3 and Oi(0) = 0. 
From Eq. (E3), we obtain 

- i cos^ kr^k^P' + [\kr^ cosfcr* sinfcr, + \ sin^ krA kP = S{k) , (28) 
9 \ 9 o J 



where r» = Csi]^ and 



S{k)=ATT i C{r)sinkrdr. (29) 







Now let us describe our method. When the data of C(r) are given, S{k) is computed from the Fourier 
sine transform of it. In actual situations, C(r) is to be obtained by observation. However, here we use 
C{r) obtained from Eq. ( [l^ ) by giving an original spectrum by hand. Then we solve Eq. ( p8| ) and compare 
the solution with the original spectrum. Note that we do not use Eq. (Km for the evaluation of C{r) 
since it is an approximate formula obtained by taking the limit d —> 00. This is because our purpose is 
to examine the accuracy of our method which uses several approximations. For the distance to the LSS, 
we choose d— 100 r*. 

The theoretical expression for C(r) includes contributions from Ci of all £; < I < 00. In reality, 
the monopole and dipole contributions of the anisotropy cannot be observed and the higher multipole 
moments are limited by the angular resolution of observation. Therefore, we subtract these contributions 
from C{r) and define the following function instead of Eq. (||), 

a,.(r) = ^^Gfl 1-^ , (30) 
1=2 ^ ^ 

where Iq is an upper limit of I and we take Zq ~ 2000 in the following calculation. We denote the 
corresponding C(r) by Cohsij')- Since the Z = 0, 1 terms contribute to C(r) mainly at large r, the source 
term is not expected to be modified significantly by neglecting these terms. As for the contributions 
of Z > Zq, they are negligible since C; decreases exponentially for large I. In Fig. |l|, we show Cir) and 
Cohs{T^ for a power law spectrum with a damping factor, P(fc) = (fcc?)"'^ exp(— fc/^o). In this case, C{r) 
can be calculated analytically and we can estimate the effect of using the finite data set of C; . We also 
plot Capp{r) which is given by cutting off the high multipole contributions of I > Iq from C(r). The 
relative error between Capp{r) and the analytic C{r) is below 2 %, which justifies the neglect of the high 
multipole contributions. On the other hand, Cobsif) differs from C{r) substantially. However, since we 
are interested in scales much smaller than d, this subtraction of / = 0, 1 modes is expected to be (and in 
fact found to be) harmless. 
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FIG. 2. The exact source term (solid line) and the source term calculated from the Fourier transform of Cobs{r) 
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(dotted line) for the spectrum same as Fig. |l|. The relative error, which results from the small scale approximation, 
r/d and the finiteness of the C'l data, is below 5%. 



For a given spectrum, the exact S{k) can be evaluated by simply substituting the original P{k) into 
the left-hand side of Eq. (28). Since we are interested in the inversion problem, we need to know how 
accurately the source term S{k) can be calculated from Cotsif)- In particular, we must reproduce the 
exact one with a good accuracy on scales I > 100. We therefore examine the accuracy of the calculated 
S{k) by comparing it with the exact one. 

Before we make this comparison, however, there is one more problem to be resolved. It is the problem 
of the range of definition of r. The Fourier transform of C(r) cannot be performed in the exact sense 
since r given by Eq. ( p^ ) is not defined for r > 2d. To avoid this difficulty, we cut off Cobs (?") by a smooth 
function which damps out steeply at some scale r ^ O.Old which is much smaller than 2d but large 
enough to reproduce the spectrum over a sufficiently wide range. Although this lowers the amplitude of 
large fc-modes, the function obtained by the Fourier transform has little deviation from the exact S{k) on 
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scales of our interest. Both the exact and calculated S{k) are shown in Fig. 0. Since we take Iq — 2000, 
the region where the calculated source term matches the exact one is below Ka — I < 2000. 

As mentioned before, Eq. (Eq) has singularities at fcr* = {n + l/2)7r which cause difficulties when we 
solve it numerically. We can, nowever, determine the values of P{k) at these singularities if we assume 
that the derivative of P{k) is finite. Then the first term on the left-hand side of Eq. (^ ) vanishes at 
fcr* — {n + ^)7r, and the values of P{k) at these singularities are given by 
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Once these values are given, we can solve Eq. ( p8| ) by expanding it around the singularities. We search the 
true solution which connect the adjacent singularities using the shooting method. We solve the equation 
until the 5th singularity, kd — 4507r, for the following two cases of the original spectra: 
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(32) 
(33) 



The first one is a double power-law spectrum and the second is a single power-law spectrum with a spiky 
structure, either with a peak (s = -1-1) or a dip (s = —1). The results for the choice of the parameters as 
p = 2, A = 10, kpd = 600, cr = 10 and kgd = 1000 are shown in Fig. ||. We find our method reproduces 
the original spectra with a good accuracy. In particular, even if the spectrum has a sharp peak or a dip, 
we can resolve such a local structure using this method. The numerical solution diverges as it approaches 
the singularities (indicated by the triangles), but the relative error except for the regions close to the 
singularities is below 5%. 

For comparison, in Fig. ^, we also show the corresponding CMB angular power spectra calculated by 
using the transfer functions given by Eq. (p7[). In the case of the spectrum with a peak, one can clearly 
see an enhancement of C; at / ~ 600, thougnthe peak height is suppressed by a factor of '--^ 4 as compared 
to the original peak in the primordial spectrum. As for the spectrum with a dip, the structure is even 
more suppressed and leaves only a small imprint in the angular spectrum. In a realistic case, a small 
structure may be present in the CMB spectrum as a consequence of an irregular feature in the primordial 
spectrum. The above result indicates that our method is capable of reproducing such a feature from the 
CMB spectrum. 
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FIG. 3. The original spectra (the solid curves) and the reconstructed spectra (the boxes and triangles). The 
left panel shows the case of the double power-law spectrum given by Eq. (pd) and the middle and the right panels 
show the single power-law spectra with a sharp peak (s = +1) and a dip (s = —1), respectively, given by Eq. (p^). 
The triangles show the locations of the singular points. We stopped the numerical integrations in the vicinity of 
the singularities when the relative error exceeded 10 %. 
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FIG. 4. The CMB angular power spectra for the two initial spectra, Eq. (^) (left panel) and Eq. (^3|) (right 
panel). The solid line and the dashed line on the right panel show the case with a peak (s = +1) and the case 
with a dip (s — —1), respectively, and the dotted line shows the case without the spiky structure in the initial 
spectrum. The normalization is arbitrary. Since we have neglected the contribution of baryons to the transfer 
functions, the acoustic peaks appear less prominent and their locations are shifted to larger I. For realistic models 
with baryons, the acoustic peaks will become more prominent and their locations will shift to smaller I. For 
example, the first peak, which is located at I ~ 300 in the figures, will shift to I ~ 200. 



Finally, it is worthwhile to comment that the presence of the singularities in the differential equation 
( P3^ ) may be regarded as an advantage, since the values of P{k) at the singularities can be estimated 
without solving the differential equation. In particular, if there is a good reason to believe that the 
spectrum should be a smoothly varying function, a qualitative feature of the spectrum can be obtained 
at once. For example, in the case of the double power-law spectrum (^2[), one can see that the original 
spectrum can be approximately recovered by simply interpolating between the adjacent triangles shown 
in Fig. I 



IV. CONCLUSION 



We have considered the problem of reconstructing the initial power spectrum of metric perturbation 
P(k) from Ci data. As a first step, we have investigated a simple case, namely, the Einstein-de Sitter 
universe with negligible baryons and negligible thickness of the LSS. In this toy model, the observed 
temperature fluctuations are represented only by the perturbation variables on the LSS and the ISW 
effect is absent. Then, the relation between the initial spectrum and the angular correlation function 
is expressed in terms of an integral equation. We have shown that this equation can be transformed 
to a second order differential equation for P{k). We have found that by forming an appropriate linear 
combination of the angular correlation function and its derivative, the order of the differential equation 
can be reduced to the first order. The resulting equation is found to have singularities that come from 
the acoustic oscillations of photons, hence their presence is inevitable in any cosmological models not 
restricted to our simple model. Fortunately, however, the presence of these singularities turns out to be 
not only harmless but rather advantageous. At the singularities the coefficient of P'(fc) vanishes, and the 
values of P{k) there can be obtained without solving the differential equation. By plotting the values 
of P{k), one can obtain a rough estimate of the behavior of the spectrum. Then, to recover the precise 
features of the spectrum, the rest part of P{k) can be solved with its values at the singularities as the 
initial values. We have found our method can reproduce the original spectrum with a good accuracy even 
for a spectrum with a sharp, spiky structure. 

The method presented here is applicable only to the Einstein-de Sitter universe in which the ISW 
effect is negligible. The ISW effect gives an important contribution even in a flat universe model if 
the cosmological constant is present. Such a model, called the ACDM model is preferred by recent 
observations |13|. Our next step is to include the ISW effect in our method which is currently under 
study. 
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